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Abstract 

We review the Laporta algorithm for the reduction of scalar integrals to the 
master integrals and the differential equations technique for their evaluation. 
We discuss the use of the basis of harmonic polylogarithms for the analytical 
expression of the results and some generalization of this basis to wider sets of 
transcendental functions. 
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1 Introduction 



In the last few years a large amount of work has been devoted to the improvement of 
the techniques for the calculation of Feynman diagrams. The reason is that future 
high energy physics experiments will reach a measurement precision that will require, 
from the theoretical counterpart, the control on the NNLO quantum corrections for 
several physical observables. 

Basically two approaches have been developed for the calculation of Feynman 
diagrams: the first one is based on the numerical and the other on the analytical 
evaluation of the integrals involved. The goal of both approaches is the complete and 
"automatic" evaluation of Feynman diagrams in multi-scale processes, but, nowa- 
days, this goal is far from being achieved. Different problems arise. While for the 
numerical approach the presence of different scales is not a problem, the treatment 
of infrared singularities, thresholds and pseudo-thresholds is of complicated solution 
and sometimes it has to be performed in a semi-analytic way. On the other hand, 
the analytic approach gives a complete control on the "difficult" regions of the spec- 
trum, but it is, at the moment, constrained to processes in which the scales in the 
game are at most three. 

Nevertheless, in both cases great results have been obtained. 

In a semi-numerical approach to the calculation of two-loop Feynman di- 
agrams was proposed and applied to the two-loop self-energy of the Higgs boson 
and to the decay Z — > bb; in |2j a method based on the Bernstein-Tkachov the- 
orem was proposed and applied both to multi-leg one-loop and to two-point and 
three-point two-loop Feynman diagrams. In [3] a numerical method based on the 
sector decomposition was proposed and applied to multi-leg one-loop Feynman dia- 
grams calculations as well as to two-loop and three- loop two-, three- and four-point 
functions in the non-physical region and to the evaluation of phase-space integrals. 
Finally, in [3] the numerical evaluation of two-point functions was made by means 
of differential equations solved with the Runge-Kutta method. 

For what concerns the analytical approach, the evaluation of vacuum diagrams, 
two-point, three-point and four-point functions (see for example [5], [51 ITU ITo] 127]. 
and Uni I2H1 ) j used also in the case of phase-space integrals (see [8j, was made, in 
the last few years, with a variety of different techniques. Since jU], the most used one 
consists in the reduction of the Feynman diagrams to a small set of scalar integrals, 
via integration-by-parts, Lorentz-invariance [TD] and general symmetry identities 
followed by the calculation of the scalar integrals with different methods such 
as, for example, the expansion by regions ^2], the Mellin-Barnes transformations 
|13j . the relations among integrals of different dimension D [14j, or the differential 
equations method [To] . 

In this paper we will review the algorithm for the reduction of the Feynman 
diagrams to the set of independent scalar integrals, called Master Integrals (Mis) 
and the differential equations technique for their evaluation. The problem of the 
choice of the basis of functions used for the expression of the analytical results will 
be also discussed, giving particular emphasis to the Harmonic Polylogarithms and 
the several extensions that took place recently. 
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Figure 1: Flowchart of the method used for the reduction to the Mis. 

2 Algebraic Reduction to Master Integrals 

The calculation of a physical observable for a certain reaction in perturbation the- 
ory is connected to the evaluation of the Feynman diagrams involved in the process. 
Once the observable is written in terms of Feynman diagrams and the Dirac algebra 
is performed (that means usually that the traces over the Dirac indices are eval- 
uated), we find an expression which is a combination of (several) scalar integrals, 
whose ultraviolet and infrared divergences are regularized within the dimensional 
regularization scheme. The general structure of such integrals at the 2-loop level is 
the following: 



T> ni . . . T) nt ' 



where t is the number of different denominators T>i (this number is called topology of 
the integral), i is the number of independent scalar products Si on the numerator, 
and D D k stands for a suitable integration measurement (normalization) for the 
integrals. 

The reduction of the integrals of Eq. to the Mis is schematically represented 
in the flowchart of Fig. ^ and it is based on the following steps: 

1. Once the Feynman diagrams for the evaluation of the observable are written, 
we project the corresponding amplitude on a basis of known tensors in such a 
way that the amplitude can be written in terms of scalar form factors. These 
form factors are expressed in terms of a huge number of scalar integrals, that 
can have on the numerator a combination of scalar products between an ex- 
ternal momentum and a loop-momentum or between two loop-momenta; note 
that, because we use dimensional regularization, a simplification between a 
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scalar product and a denominator that contains this scalar product is possi- 
ble. Once a denominator disappears, because of the simplification against a 
scalar product, the topology of the integral is of course lowered by a unit. In 
principle, given an integral of topology t, one must consider all the possible 
t\ subtopologies found simplifying repeatedly a denominator against a scalar 
product on the numerator. 

2. The scalar integrals are classified with respect to their topology. Because we 
are dealing with integrals that have a suitable mass distribution (for exam- 
ple we can consider integrals with mass-less propagators and outgoing legs) it 
can happen that subtopologies coming from different simplifications are the 
same subtopology, but expressed with different routings. Therefore, the "inde- 
pendent" subtopologies have to be chosen and the "dependent" subtopologies 
have to be transformed by means of suitable transformations in the indepen- 
dent ones. This analysis can be done in the framework of the "Auxiliary 
diagram scheme" or the "shifts scheme", as explained extensively in j!6j . 

3. The reduction of the scalar integrals belonging to the independent topologies 
to a "hopefully" small set of master integrals (Mis), is done by means of the 
so-called Laporta algorithm ^7j, which consists in the following. We know 
that the .D-regularized scalar integrals coming from the projection operation 
are not all independent, but they satisfy certain classes of identity relations. 

• The most important class is constituted by the integration-by-part iden- 
tities (IBPs), introduced in IBPs link scalar integrals of the same 
topology, but with different power of the denominator and different scalar 
products on the numerator, among each other and to scalar integrals of 
subtopologies. IBPs can be written in the following way: 

where i = 1, 2 and where t> M = k% , k%,pf is one of the independent vector 
of the problem. In the case of a 3-point functions, Eq. (j2J) gives 8 equa- 
tions for initial scalar amplitude in the brackets. For a 4-point function, 
instead, the IBPs are 10. 

• Another class of identity relations that can be used in the reduction pro- 
cess is related to the fact that the integrals are Lorentz scalars ^Hj ■ This 
property translates into one additional equation for a 3-point function: 

2 

(piP2 - v\vi) 

n=l 

where p± and P2 are the independent vectors of the problem, or three 
additional equations for a 4-point function: 
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(P2P3 ~ P2P3) Yl 
n=l 

where the independent vectors are pi, P2 and p%. 

• In the case in which a suitable mass distribution is considered, we can 
find additional equations considering the symmetries of the diagrams (see 
jTTj). In general a symmetry of the problem brings to an identity that 
the integrals have to satisfy. 

Considering all the identity relations mentioned above, a huge system of linear 
equations is constructed. The unknowns are the scalar integrals themselves 
and the point is that the construction of the system can give more equations 
than unknown amplitudes ^7j, overconstraining formally the system; but not 
all the equations are independent. It can happen that the system is effectively 
overconstrained, and then the solution of the system gives the integrals of the 
topology under consideration as a combination of the Mis of the subtopologies, 
or it is not, and then all the integrals of the topology under consideration are 
expressed as a combination of the Mis of that topology (and the Mis of the 
subtopologies). The solution of the system is performed with the Gauss law of 
substitution. The entire chain is completely algebraic and can be implemented 
in a computer program. 

Several authors developed own programs, written in FORM [18 , C, or Mathematica 
[TT?] , for the generation of the linear system and for its solution EE| • Recently, 
a computer program written for Maple [22], using the Laporta algorithm for the 
reduction of scalar integrals to the Mis was published [23]. 

3 The Differential Equations Method 

The Mis are functions of the external kinematical invariants and, therefore, they 
satisfy a system of first-order linear differential equations. As a simple example, 
consider the case in which a topology t has a single MI. Let us choose the basic 
scalar integral of the topology 3 : 

F( Si ) = J D%D% — 1 — , (7) 

3 The choice of the set of Mis to which all the scalar integrals belonging to a certain topology can 
be reduced is totally free. One choice is connected to another by the identity relations constructed 
for the reduction process. A criterion in the preference of one set with respect to the others is, 
of course, related to the solution of the system of differential equations. We choose the set that 
satisfy the easiest possible system. That means, usually, the set for which the system is decoupled 
exactly in D or at least triangularizes in the limit D — > 4. 
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where Sj stand for the independent invariants that can be constructed with the 
external momenta of the problem (for example s\ = p\, S2 = pi ■ P2, etc.). 
The following matrix can be constructed: 



As F(si) depends on the invariants Sj, we have, on one hand: 

OjM £-F( ai ) = J2 a U k (si)^F( Sl ) } (9) 

where the functions a£ jy(s») are linear combinations of Sj. On the other hand, 
because F(si) is evaluated in dimensional regularization, and then it is convergent, 
we can perform the derivative of Eq. (jHJ) directly on the integrand ../V it1 
getting a combination of scalar integrals with additional power on the denominator 
and scalar products on the numerator. Using the identity relations loaded for the 
reduction to the Mis, the system of Eq. (jHJ) becomes a linear system involving the 
derivatives of F(si) with respect to the invariants Sj, F(si) itself and the Mis of the 
subtopologies. The system can be inverted and the differential equations can be 
written in the following general way: 

■^rF(si) = A(D, Si ) F( 8i ) + n(D, Si ) , (10) 

where the non-homogeneous term Q(D, Si) contains the Mis of the subtopologies, 
that have to be considered known. The homogeneous term A(D, s;) gives the analyt- 
ical structure of the function F(si), containing the thresholds and pseudo-thresholds 
of the diagram under consideration, that appear as singular points for the differential 
equation 4 . Note that one of the differential equations is sufficient for the solution of 
F(si), provided that we are able to fix the boundary conditions. Therefore, we can 
consider only the equation with respect to, say, s. The solution of Eq. (jTUJl is done 
by means of the Euler's method of variation of the arbitrary constants. 
We can sketch the search for the solution as follows: 

1. We expand F(s) (the dependence on the other invariants is understood) in 
Laurent series of (D — 4): 

n 

F(s) = J2 (D - 4)^( a ) , (11) 

j=-k 

where k is the maximum pole and n is the required order in (D — 4) needed 
for F(s). Order by order in (D — 4) we have, then, to solve the equation: 

^F j (s) = A(s)F j (s)+n j (s), (12) 



4 The construction of the system of first-order linear differential equations is outlined in the 
flowchart of Fig. 



5 



Master 


Integrals 






r- 


- Generation of Diff. Eqs. (+ IBP's) 


System 


of Diff. Eqs. 





Solution of the System in Laurent series of e 



Figure 2: Flowchart of the Differential Equations method. 

where for j = —k (the maximum pole) we have f2_^(s) = fi_^(s) (f2(s) is the 
non- homogeneous part in Eq. ([10))) and for j > —k, Cl_k(s) can involve also 
the previous orders in (D — 4) of F(s). Note that, while the non-homogeneous 
term fij(s) is different at each order in (D — 4), the homogeneous equation is 
the same. 



2. We solve the homogeneous equation 



§- s m = A( S )f(s), 



finding the solution: 



f(s) = exp 



A(t)dt 



(13) 



(14) 



3. We express the solution of the non- homogeneous equation order by order in 
(D — 4) in the integral form: 



F 3 {s) 



i 



/(*) 



Clj(t)dt + kj 



(15) 



where kj is the arbitrary constant of integration. 



4. We fix the constant of integration imposing the initial conditions. In order 
to find the initial conditions we have to know additional pieces of information 
about the integral we are calculating. For example it is sufficient to know that 
the integral is regular for some value of s. 

In the case of N Mis (N > 1) we have still a system of iV first-order coupled linear 
differential equations for every variable Sj. As in the previous case the equations in 
one of the invariants are sufficient for the solution. 

In spite of the elegance of the method, two problems arise. 

One is connected with the number of Mis that a topology can have. In fact, 
while the solution of a first-order linear differential equation (case with one MI) 
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is relatively trivial, a second-order differential equation (case with two Mis) can 
give more problems and starting from the third-order one it can be hard to find 
the solution, except in very particular cases. We can understand, therefore, the 
importance of the choice of the set of Mis. A choice that, even in the case of two or 
more Mis, could triangularize the system, at least in the (D — 4) expansion, would 
be of course more suitable than another; but this choice can not always be done 
and, in general, there is not always a solution for the problem of a topology with 
many Mis. 

Another problem concerns the choice of the basis of functions used for the ex- 
pression of the analytical results. Even in this case, there is not a unique solution, 
but, nevertheless, one can follow some guidelines. For example the uniqueness of 
the representation of the result in terms of these functions; the non-redundancy of 
the representation and the absence of "hidden zeroes" (that means to avoid repre- 
sentations that can give expressions that are not manifestly zero, but, because of 
the fact that these functions satisfy certain relations, they are effectively zero); the 
total control on the expansions in all the points of the domain and on the analytical 
continuation (these two "properties" are needed for a numerical evaluation of the 
functions and then of the result). 

All these properties are fulfilled, by construction, by the set of functions called 
Harmonic Polylogarithms, that we are going to review briefly in the next paragraph. 

4 Harmonic Polylogarithms and Related Gener- 
alized Functions 

In many cases, a suitable integral representation for the Feynman diagrams can be 
found in terms of hypergeometric functions. Nevertheless, in problems involving 
many scales, the dependence of the generalized hypergeometric functions on these 
scales is highly non-trivial. Moreover, if we regularize the divergent integrals in 
dimensional regularization, an expansion in (D — 4) is required, for renormalization 
purposes. The expansion of a hypergeometric function in its parameters can be very 
complicated j21]. This is the reason why one looks for a solution of the differential 
equations directly expanded in (D — 4). Order by order in (D — 4), we can find a 
representation for the Mis in terms of Nielsen's polylogarithms and related functions 
|25j . But, firstly, the representation in terms of polylogarithms suffers from the 
problem of "hidden zeroes" discussed at the end of the previous paragraph; moreover, 
it can happen that one needs an integral representation that does not belong to this 
class of functions. 

One of the most elegant solutions for these two problems is the introduction of a 
basis of transcendental functions defined by repeated integration over a set of basic 
simple functions. Let us consider for example a case in which we have at most two 
different scales in the problem, say s and a, in such a way that we can construct 
the dimensionless parameter x = s/a; moreover, the structure of the thresholds of 
the Feynman diagrams involved in the calculation is such that the possible singular 
points are x = and x — 1 and no squared roots are present in the homogeneous 
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equation for the Mis (this is actually the case of the majority of the calculations in 
HU Uni 123 )• In this situation a suitable basis of functions in which the results 
can be expressed are the 1-dimensional harmonic polylogarithms (HPLs) introduced 
in We consider the following set of three basis functions: 

g(0;x) = -, g(l;x) = —^—, g(-l; x) = — |— , (16) 

x 1 — X 1 + X 

we define the HPLs of weight 1 as: 

r x dt 

H(0,x) = J y = logx, (17) 

r dt 

H(l,x) = I T — = -] og (l-x), (18) 

H(-l,x) = J — = log(l + a;), (19) 

and we iterate the previous definition introducing the HPLs of weight w + 1 

H(0 W+1 ; x) = ] \og w+1 x , H(a, w- x) = [ /(a; x')H(w, x'), (20) 

where a can take the values -1, and 1 and w is a vector with w components, 
consisting of a sequence of -1, and 1 as well. 

The set of functions so defined satisfy two important properties: i. they form 
a shuffle algebra, in which the product of two HPLs of weight W\ and W2 is a HPL 
of weight Wi + W2', ii- they form a closed set under certain transformations of the 
argument. In particular, this property is needed for the asymptotic evaluation of 
such a functions, connecting the expansion in x — > to the one in x —>■ 00. Another 
important property of HPLs is that their analytical structure is manifestly shown. 
The integral representations of Eq. (fTK|) can be expressed in terms of HPLs directly 
or re-conducted to HPLs via simple integration by parts. 

The choice of the basis f|16|) for the expression of the HPLs is connected with 
the structure of the thresholds of the diagrams involved in the calculation. It can 
happen that the solution of the homogeneous differential equation contains a square 
root, that can not be avoided by a suitable change of variables. Moreover, we can 
deal with a calculation in which, for example, the diagrams can have three different 
kind of thresholds: x — 0, x — 1 and x = 4. In this case an extension of the basis of 
HPLs is needed j^ZJ- We introduce the following additional basis functions: 

3(:p4;x) = — — , g(c,x) = p ?f , g(c,x) = -= 

x 2 2 2 ~ 2 

1 _ 1 

g(^r,x) = , #(=Fl=Fr,x) 



y/x{4 ± X) ' y/x(4 ±X) (1 ± x) 

g (±l/A;x) = T ^—, g(±l±r/4;x) ' 
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g(r /A;x) = i , g(-r /A;x) = 1 ■ (21) 

The relative weight- 1 HPLs are defined as the integration between and x of 
the previous basis functions (all the functions in Eq. (|21|) are integrable in x = 0) 
and the generalization to higher weights is done by Eq. ([2*0]) in which now the single 
weight can take the values 0, ±1, ±4, ±r, c, c, =Fl =f r, ±1/4, ±1 ± r/4, ±r /4. 

The set of generalized HPLs so introduced satisfy all the properties of the former 
set of HPLs. Similar generalizations can be carried out in the case of different 
thresholds. 

Another generalization can be done in the case in which we have three scales 
in the game and, therefore, two dimensionless parameters can be constructed. This 
gives the so-called 2-dimensional HPLs introduced in [2B1 [29j and their further ex- 
tension introduced in [30] . 

As a remark, note that the representation of the analytical results in terms of 
HPLs allows a perfect control on the numerical evaluations. 

5 Summary 

In this paper the Laporta algorithm for the reduction of the Feynman diagrams 
to the master integrals and the differential equations technique for their analytical 
evaluation are reviewed. A particular attention is paid to the problem of the choice 
of the basis of functions in which the result can be expressed. Some argument for 
the use of HPLs is done. In particular, the possibility of an extension of the set of 
functions to the case of problems involving different thresholds and different scales 
(more than two) is briefly outlined. 
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